# ==============================================================================
# DESCRIPTION ----
# ==============================================================================

# Conducts differences-in-differences analysis of effect of former 
# bureaucrats on NPO's negotiated contract values,
# using yearly aggregated data. 

# ==============================================================================
# PREAMBLE, LIBRARIES, AND IMPORT ----
# ==============================================================================
# Preamble, packages -----------------------------------------------------------
library(tidyverse)
library(DIDmultiplegt)
library(fect)
library(panelView)
library(benford.analysis)
library(gridExtra)
library(plm)
library(tsibble)

# Import user-created functions
source("code/0.2_functions.R")

# Import subsidy data
load("data/jnpo.Rdata")

# ==============================================================================
# TRANSFORM SUBSIDY DATA TO YEARLY TIME SERIES ----
# ==============================================================================

# Transform into monthly time series dataset -----------------------------------
# Collapse data to firm year-month level: without ministry information
npo_amakudari <- jnpo %>%
  filter(!is.na(govt_reemployees), amount > 0, grant_year > 2007) %>%
  mutate(amakudari = as.numeric(as.character(govt_reemployees))) %>%
  group_by(grantee_clean, grant_year) %>%
  summarize(amount = sum(amount), amakudari_n = max(amakudari)) %>%
  mutate(
    amakudari = ifelse(amakudari_n > 0, 1, 0),
    amount_log = log(amount),
    amount = amount / 1000000
  )

# ==============================================================================
# ANALYSIS ----
# ==============================================================================
#### FEct, IFECT, and MC methods -----------------------------------------------
# FEct
set.seed(123)
fect <- fect(amount_log ~ amakudari, data = npo_amakudari, 
             index = c("grantee_clean", "grant_year"), 
             method = "fe", force = "two-way", na.rm = TRUE, se = TRUE)

fect_plot <- plot(fect, main = "FEct", 
                  ylab = "Effect on log(contract value)",
                  cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, proportion = 0.25)

set.seed(123)
fect_level <- fect(amount ~ amakudari, data = npo_amakudari, 
                   index = c("grantee_clean", "grant_year"), 
                   method = "fe", force = "two-way", na.rm = TRUE, se = TRUE)

fect_plot_level <- plot(fect_level, main = "FEct", 
                        ylab = "Effect on contract value", ylim = c(-500, 500),
                        cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, proportion = 0.25)

# IFEct
set.seed(123)
ifect <- fect(amount_log ~ amakudari, data = npo_amakudari, 
              index = c("grantee_clean", "grant_year"), 
              method = "ife", force = "two-way", na.rm = T, se = TRUE)

ifect_plot <- plot(ifect, main = "IFEct", 
                   ylab = "Effect on log(contract value)", 
                   cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, proportion = 0.25)

set.seed(123)
ifect_level <- fect(amount ~ amakudari, data = npo_amakudari, 
                    index = c("grantee_clean", "grant_year"), 
                    method = "ife", force = "two-way", na.rm = T, se = TRUE)

ifect_plot_level <- plot(ifect_level, main = "IFEct", 
                         ylab = "Effect on contract value", ylim = c(-500, 500),
                         cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, proportion = 0.25)

# Matrix completion
set.seed(123)
mc <- fect(amount_log ~ amakudari, data = npo_amakudari,
           index = c("grantee_clean", "grant_year"), 
           method = "mc", force = "two-way", na.rm = T, se = TRUE)

mc_plot <- plot(mc, main = "Matrix Completion", 
                ylab = "Effect on log(contract value)",
                cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, proportion = 0.25)

set.seed(123)
mc_level <- fect(amount ~ amakudari, data = npo_amakudari,
                 index = c("grantee_clean", "grant_year"), 
                 method = "mc", force = "two-way", na.rm = T, se = TRUE)

mc_plot_level <- plot(mc_level, main = "Matrix Completion", 
                      ylab = "Effect on contract value", ylim = c(-500, 500),
                      cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, proportion = 0.25)

plot_binary <- arrangeGrob(
  fect_plot, ifect_plot, mc_plot, fect_plot_level, ifect_plot_level, mc_plot_level,
  ncol = 3,
  top = "Estimated ATT by estimator and log or level (in million yen) outcome, with binary treatment"
)

ggsave("figures/a35_npo_robust_yearly.pdf", plot_binary, height = 6, width = 10)

# Classic TWFE -----------------------------------------------------------------
# Binary amakudari indicator (log outcome)
fit_binary <- plm(amount_log ~ amakudari,
                  data = npo_amakudari, 
                  index = c("grantee_clean", "grant_year"), 
                  model = "within", 
                  effect = "twoways")

# Binary amakudari indicator (level outcome)
fit_binary_level <- plm(amount ~ amakudari,
                        data = npo_amakudari, 
                        index = c("grantee_clean", "grant_year"), 
                        model = "within", 
                        effect = "twoways")

# Continuous amakudari indicator (log outcome)
fit_continuous <- plm(amount_log ~ amakudari_n,
                      data = npo_amakudari, 
                      index = c("grantee_clean", "grant_year"), 
                      model = "within", 
                      effect = "twoways")

# Continuous amakudari indicator (level outcome)
fit_continuous_level <- plm(amount ~ amakudari_n,
                            data = npo_amakudari, 
                            index = c("grantee_clean", "grant_year"), 
                            model = "within", 
                            effect = "twoways")

# Create table for classic 2WFE
stargazer::stargazer(fit_binary, fit_continuous, fit_binary_level, fit_continuous_level,
                     omit = "grantee_clean", 
                     omit.stat = c("adj.rsq", "rsq"),
                     type="latex",
                     out = "tables/a23_twfe_yearly.tex")

#### DiDm ----------------------------------------------------------------------
# Define parameters
Y = "amount_log"
G = "grantee_clean"
T = "grant_year"
D_binary = "amakudari"
D_cont = "amakudari_n"

# Run and plot models
set.seed(123)
par(mar=c(4,4,4,4), mfrow=c(2,1), cex.axis=0.5, cex.main=0.75)
neg_binary <- did_multiplegt(npo_amakudari, Y, G, T, D_binary, placebo = 3, brep = 100)
abline(h=0, col="black", lty = 2)
title(main = "Bureaucratic hires (binary)")
neg_cont <- did_multiplegt(npo_amakudari, Y, G, T, D_cont, placebo = 3, brep = 100)
abline(h=0, col="black", lty = 2)
title(main = "Bureaucratic hires (continuous)")

# FIGURE A34
# Plot models: binary
binary_plot <- as_tibble(neg_binary) %>%
  select(-contains("N"), 
         effect_1 = placebo_1, effect_2 = placebo_2, effect_3 = placebo_3,
         se_effect_1 = se_placebo_1, se_effect_2 = se_placebo_2, 
         se_effect_3 = se_placebo_3) %>%
  pivot_longer(cols = c(starts_with("effect"), starts_with("se"))) %>%
  mutate(
    effect = ifelse(str_detect(name, "se"), "se", "effect"),
    time = c(-3, -2, -1, 0, -3, -2, -1, 0)) %>%
  select(-name) %>%
  pivot_wider(names_from = effect, values_from = value) %>%
  mutate(model = "Bureaucratic hires (binary)")

# Continuous plot
cont_plot <- as_tibble(neg_cont) %>%
  select(-contains("N"), 
         effect_1 = placebo_1, effect_2 = placebo_2, effect_3 = placebo_3,
         se_effect_1 = se_placebo_1, se_effect_2 = se_placebo_2, 
         se_effect_3 = se_placebo_3) %>%
  pivot_longer(cols = c(starts_with("effect"), starts_with("se"))) %>%
  mutate(
    effect = ifelse(str_detect(name, "se"), "se", "effect"),
    time = c(-3, -2, -1, 0, -3, -2, -1, 0)) %>%
  select(-name) %>%
  pivot_wider(names_from = effect, values_from = value) %>%
  mutate(model = "Bureaucratic hires (continuous)")

plot <- rbind(binary_plot, cont_plot)

# Create combined plot
ggplot(plot, aes(time, effect)) +
  geom_point(color = "steelblue2", size = 1.5) + 
  geom_line(color = "grey80") +
  geom_errorbar(aes(x = time, 
                    ymin = effect - 1.96*se, 
                    ymax = effect + 1.96*se),
                color="grey30", size=0.5, alpha = 0.5, width = 0.1) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  ylab("Effect of appointment on log(contract value)") + 
  xlab("Time (periods before hire)") +
  scale_y_continuous(limits = c(-1.5, 1.5), breaks=seq(-1.5, 1.5, 0.5)) +
  facet_wrap(~model) +
  theme_classic()

ggsave("figures/a34_npo_amakudari_yearly.pdf", height = 3.5, width = 7)
